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Preface 


1  his  report  describes  research  into  several  areas  of  groundwater  model¬ 
ing.  The  primary  focus  of  this  work  is  basic  research  on  computational  al¬ 
gorithms,  but  a  secondary  aspect  is  a  study  of  the  use  of  scientific 
visualization  technology.  The  work  was  funded  unde:  the  In-House  Labo¬ 
ratory  Independent  Research  (ILIR)  program  at  the  U.S.  Army  Engineer 
Waterways  Experiment  Station  (WES).  Vicksburg,  MS.  Additional  re¬ 
search  and  the  publication  of  this  report  were  funded  by  the  Strategic  Envi¬ 
ronmental  Research  and  Development  Program  (SERDP),  Department  of 
Defense. 


Dr.  James  H.  May  and  CPT  Edward  Mazion,  Engineering  Geology 
Branch,  Eartho,uake  Engineering  and  Geophysics  Division,  Geotechnical 
Laboratory  (GL),  WES,  performed  the  2-D  Rocky  Mountain  Arsenal 
(P.MA)  study  described  in  thi.s  report.  Mr.  John  B.  Palmerton,  Rock  Me¬ 
chanics  Branch,  Soil  and  Rock  Mechanics  Division,  GL,  WES,  performed 
the  analysis  of  Cerrillos  Dam  in  Puerto  Rico  discussed  in  this  report. 

Most  of  the  scientific  visualization  work  for  RMA  was  done  by  Mr.  Scott 
Weberg,  Scientific  Visualization  Center  (SVC),  Advanced  Technology 
Center  (ATC),  Information  Management  Division  (IMD),  Information 
Technology  Laboratory  (ITL),  WES.  Mr.  Charles  S.  Jones,  SVC,  ATC, 
IMD,  ITL,  WES,  assisted  in  the  visualization  for  contaminant  transport. 
Cray  Y-MP  improvements  were  made  by  Mr.  Alex  Carrillo,  Department  of 
Defense  High  Performance  Computing  Center,  ATC,  IMD,  ITL,  WES. 

Dr.  Fred  T.  Tracy,  Interdisciplinary  Research  Group,  Computer-Aided  En¬ 
gineering  Division  (CAED),  ITL.  WES,  wrote  this  report  (Mr.  Carrillo  pro¬ 
vided  an  excellent  draft  documenting  his  work)  and  performed  the 


remaining  work  described  herein.  Dr.  iNccu  l^.  iviusnci  ^^as  r\ctiiig  v^itici 


of  CAED  during  this  study  and  preparation  of  this  report,  and  Director  of 


ITL  was  Dr.  N,  Radhakrishnan. 


This  work  was  coordinated  with  the  WES  Groundwater  Modeling  Team 
who.se  chairman  is  Dr.  Jeffery  P.  Holland.  Director.  Computational  Hydrau¬ 
lics  Institute,  Hydraulics  Laboratory,  and  Program  Manager,  Groundwater 
Modeling  Program. 
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At  ihe  time  of  publication  of  this  report.  Director  of  WES  was 
Dr.  Robert  W.  Whalin.  Commander  was  COL  Bruce  K.  Howard,  EN. 


The  contents  of  this  report  are  not  to  be  used  foi  advertising,  publication, 
or  promotional  purposes.  Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
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Introduction 


Why  This  Research? 

The  flow  of  groundwater  is  an  extremely  complex  process  to  model  be  ¬ 
cause  of  the  diverse  types  of  flow  and  the  heterogeneous  nature  of  the  po¬ 
rous  iTiedia.  When  the  transport  of  diverse  types  of  contaminants  is  added 
to  the  problem,  the  difficulty  and  complexity  of  modeling  real-world  three- 
dimensional  (3-D)  groundwater  flow  increases  by  orders  of  magnitude. 
Also,  it  is  very  difficult  to  visualize  groundwater  input  and  output  into  nu¬ 
merical  programs  because  of  the  complexity  of  the  geometry  and  heteroge¬ 
neity  of  the  porous  medium.  This  report  documents  basic  and  applied 
research  in  the  area  of  groundwater  modeling  and  visualization. 


Cray  Y-MP 

The  use  of  a  supercomputer  such  as  a  Cray  Y-MP  is  also  essential  for 
the  timely  computation  of  results  of  a  large  3-D  groundwater  model.  This 
is  because  the  highly  nonlinear  nature  of  the  soil  properties  in  the  unsatu¬ 
rated  zone  creates  a  tremendous  computational  load,  and  traditional  tech¬ 
niques  do  several  iterations  at  each  time-step  in  an  attempt  to  achieve 
suitable  answers.  Further  complications  arise  because  groundwater  flow 
is  sometimes  modeled  over  extremely  long  periods  of  time,  putting  an  ad¬ 
ditional  requirement  that  the  numerics  remain  stable  and  robust.  This  re¬ 
port  describes  research  findings  in  this  area  as  well. 


Scientific  Visualization 

The  use  of  graphical  tools  to  visualize  both  the  input  data  (grid  with  ini¬ 
tial  and  boundary  conditions)  and  the  results  (pressure,  head,  velocity,  con¬ 
centration,  etc.)  is  absolutely  essential  for  successful  modeling.  This 
report  also  describes  the  scientific  visualization  technique^  developed  in 
this  research. 


Chapter  1  Introduction 
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Chapter  Summaries 


Chapter  2  uses  the  application  of  the  real-world  problem  of  flow  near 
Rocky  Mountain  Arsenal  (RMA)  to  help  answer  the  question  of  w'hen  a 
3-D  analysis  is  needed.  A  comparison  of  a  two-dimensional  (2-D)  plan 
view  time-dependent  solution  taken  to  steady  state  to  a  true  3-D  saturated 
flow  solution  where  the  phreatic  surface  is  iterated  to  the  steady- state  solu¬ 
tion  's  presented.  Both  these  solutions  use  the  finite  element  method 
(FEM).  The  details  of  what  was  done  to  create  the  3-D  FEM  program  and 
data  are  given.  Examples  of  how  visualization  was  used  are  also  given. 

The  running  times  for  the  3-D  saturated  flow  steady-slate  RMA  prob¬ 
lem  using  classical  FEM  solution  techniques  showed  that  for  3-D,  time- 
dependent,  nonlinear  problems  involving  unsaturated  or  multiphase  flow 
with  contaminant  transport,  the  running  times  even  on  the  Cray  Y-MP 
could  possibly  become  prohibitive.  Therefore,  Chapters  3  and  4  describe 
basic  research  into  otlicr  techniques  that  arc  potentially  faster  and,  for 
those  problems  where  they  can  be  applied,  should  deliver  an  equivalent 
quality  answer  much  more  efficiently.  Specifically,  the  finite  volume 
(FV)  method  is  first  applied  to  saturated/unsaturated  How  and  then  to  con¬ 
taminant  transport  and  compared  to  classical  FEM  solutions.  The  tech¬ 
niques  of  nonlinear  iteration,  operator  splitting,  etc.,  are  also  discussed. 
Both  experimental  results  and  analytical  solutions  are  used  to  compare  dif- 
ferent  procedures.  Finally,  contaminant  transport  results  are  viewed  using 
scientific  visualization  tools. 


2 


ChaptRi  1  Introduction 


2  Saturated  Flow 


Introduction 

Engineers  and  scientists  at  the  Geotechnical  Laboratory  (GL),  U.S. 
Army  Engineer  Waterways  Experiment  Station  (WES),  have  made  a  study 
of  flow  in  the  RMA  region.  The  program  used  for  the  study  is  a  modified 
version  of  a  2-D  plan  view  saturated  and  confined  flow  FEM  code 
(Warner  1987)  that  accepts  a  triangular  grid.  This  report  describes  three 
additional  things  that  were  done  as  follows: 

a.  Three-dimensional  study.  A  3-D  grid  was  built  from  the  given  2-D 
data,  and  a  3-D  FEM  unconfined,  saturated  flow  program  was  ap¬ 
plied  to  the  same  problem  with  the  2-D  and  3-D  results  then  com¬ 
pared.  This  was  done  to  help  answer  the  question  of  when  a  3-D 
solution  is  required  for  a  given  problem. 

b.  Cray  Y-MP  optimization.  Special  vectorization  techniques  were  in¬ 
vestigated  to  see  what  improvements  in  running  times  could  be 
achieved  over  the  generic  Unix  version  of  the  3-D  code.  Computa¬ 
tional  improvements  valid  for  any  scientific  computer  were  also 
made. 

c.  Scientific  visualization.  Scientific  visualization  techniques  were  ap¬ 
plied  to  the  RMA  data. 

The  3-D  computations  were  done  on  a  Cray  Y-MP  using  a  modified  ver¬ 
sion  of  a  3-D  seepage/groundwater  FEM  program  (Tracy  1991).  Scien¬ 
tific  visualization  tools  used  were  commercial  programs  such  as 
Multi-Eurpose  Graphics  System  (MPGS)  (Cray  1990)  and  customized  soft¬ 
ware  developed  in  the  Scientific  Visualization  Center  (SVC),  WES.  All 
scientific  visualization  was  done  on  Silicon  Graphics  workstations. 
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Description  of  Problem 


The  Di  oblem  consists  of  partially  confined  and  partially  unconfined 
groundwater  flow  in  a  region  near  RMA.  The  2-D  triangular  mesh  used 
for  the  pioblem  is  shown  in  Figure  1.  The  nodes  with  triangles  have 
heads  specified  as  a  boundary  condition,  and  the  nodes  identified  with  cir¬ 
cles  are  observation  wells  where  differing  amoun's  of  water  are  being  ex¬ 
tracted.  The  immediate  purpose  of  the  compu.  ton  is  a  calibration  where 
the  unknown  hydraulic  conductivities  are  adjusted  so  that  the  initial  heads 
remain  the  final  heads  after  the  problem  is  run  to  a  steady-state  solution. 

A  slurry  trench  has  been  installed  in  the  flow  region  to  modify  flow,  and  a 
zoom  of  the  grid  for  this  region  is  shown  in  Figure  2.  The  slurry  trench  is 


bEM  grid  with  boundary  conditions 


Cliapisr  2  Saturated  Flow 


Figure  2.  Zoomed  area  showing  slurry  trench 

identified  by  all  the  nodes  along  line  segment  AB.  An  impervious  wall 
such  as  a  slurry  trench  is  modeled  by  having  d.fferent  nodes  on  one  side 
of  the  wall  as  compared  to  the  other  side.  This  particular  slurry  trench  has 
many  more  nodes  on  the  left  side  than  on  the  right  side. 

There  are  two  basic  layers  which  are  alluvium  and  bedrock.  The  hy¬ 
draulic  conductivities,  however,  vary  greatly  within  these  two  broad  cate¬ 
gories  of  material.  Figure  3  shows  a  color  contour  plot  of  the  hydraulic 
conductivity  of  the  alluvium,  and  Figure  4  shows  a  color  contour  plot  of 
the  hydraulic  conductivity  of  the  bedrock.  Figure  3  shows  a  color  contour 
plot  of  the  initial  toi  il  head  values.  The  problem  is  now  completely 
specified. 
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Figure  5.  Initial  head  values 


2-D  Results 


Figure  6  shows  the  results  of  one  of  the  interinediale  calibrations.  The 
amount  of  drawdown  is  small  in  most  places,  but  there  arc  still  some  "hot 
spots”  to  fix.  It  is  clear  that  the  color  contour  plot  is  an  excellent  tool  to 
quickly  determine  the  areas  that  need  more  work.  The  final  calibration 
will  be  presented  in  the  GL  report  to  KM  A.  The  purposes  of  this  chapter 
arc  well  served  by  using  these  intermediate  results,  and  the  problem  of 
using  possibly  proprietary  information  is  avoided. 
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Figure  6.  Intermediate  calibration  drawdown  results 


3-u  Study 


At  this  point  a  study  was  begun  to  determine  how  important  a  true  3-D 
solution  is  and  how  to  visualize  the  3-D  results.  Cray  optimization  was 
also  investigated  by  personnel  in  the  DoD  High  Performance  Computing 
(HPC)  Center  group  in  ITL, 

3-D  grid 

A  3-D  grid  was  generated  from  the  2-D  data  by  adding  intermediate  lay¬ 
ers  (sec  the  hidden  line  view  of  the  grid  shown  in  Figure  7)  by  converting 
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2-D  triangles  converted  to  3-D  prisms 


Figure  9.  Specified  head  in  alluvium 
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Since  this  is  u  calibrutiun  at  a  rather  late  stage,  step  a  was  not  repeated  at 
each  iteration  in  the  3-D  groundwater  program  to  find  the  free  surface. 
However,  this  should  be  incorporated  into  the  program  for  doing  studies 
beyond  the  calibration. 

For  an  unconfined  flow  situation  where  the  water  table  is  in  the  bed¬ 
rock,  the  elevation  of  generated  .3-D  node  2  is  given  the  value  of  the  water 
elevation,  and  the  given  Q  is  divided  equally  between  generated  .3-D 
nodes  1  and  2. 


3-D  groundwater  model 

The  .3-D  groundwater  model  used  in  this  study  is  a  modification  of  one 
written  by  this  author  for  3-D  steady-state,  unconfined  saturated  flow. 

The  modifications  required  arc  us  follows; 

a.  Hydraulic  conductivity  data.  Typically,  in  an  FEM  program,  the 
material  properties  arc  specified  separately,  and  each  set  is  given  a 
number.  For  instance,  material  type  1  can  represent  sand,  and  mate¬ 
rial  type  2  can  represent  clay.  Thus,  each  element  has  a  material 
type  number  specified,  and  those  values  specified  remain  constant 
inside  the  given  element.  The  data  for  RMA  is  different.  There 
were  provided  one  hydraulic  conductivity  for  bedrock  and  one  hy¬ 
draulic  conductivity  for  alluvium  for  each  2-D  node.  Therefore,  no 
data  were  provided  for  element  properties,  but  the  node  data  were 
extended  to  include  the  two  hydraulic  conductivities.  However, 
each  element  still  had  to  be  identified  as  being  either  bedrock  or 
alluvium. 

/;.  Stiffness  matrix  computation.  Since  hydraulic  conductivity  is  not 
constant  inside  each  element,  the  stiffness  matrix  computation  had 
to  be  modified.  The  material  type  number  for  each  element  now 
represents  either  1  for  bedrock  or  2  for  alluvium.  The  finite  ele¬ 
ment  type  u.scd  in  the  3-D  model  requires  numerical  integration  to 
form  the  stiffness  matrix.  So  at  integration  point],  the  hydraulic 
conductivity  k’  is  computed  using  the  eight-node  isoparametric  ele¬ 
ment  formulation 


8 

u  =  2  v,  (1) 

/■  =  1 


The  interpolation  functions  arc 
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where 


-1  <  ^  r 

-\  <.  r\  ^  I 

-1  <  C  £  1 

\  J 

(^/i  C/)  =  n*  0  coordinates  at  node  i 

(4y,  Tjy,  t^j)  =  (^,  Ti,  coordinates  at  interpolation  point  j 

kj  -  hydraulic  conductivity  at  node  i 

c.  MFCS  files.  MPGS  files  were  written  for  the  original  grid,  the  mod¬ 
ified  grid  that  conforms  to  the  free  surface,  a  scaler  file  containing 
total  head,  and  a  vector  file  containing  Durcian  velocities. 


Visualization  of  3-0  results 

After  a  valid  3-D  grid  was  produced  and  a  correct  3-D  groundwater 
model  was  completed,  the  3-D  results  were  obtained.  It  is,  ofcour.se,  now 
very  important  to  visualize  the  resulting  voluminous  set  of  output  data.  It 
is  extremely  difficult  to  visualize  3-D  groundwater  flow  for  the  following 
three  reasons: 

u.  Relatively  flat  grid.  The  3-D  grid  in  the  /.  direction  is  extremely  flat 
compared  to  the  plan  view. 

h.  (Vo  specific  object.  In  computational  fluid  dynamics  (CFD)  there 
are  airplanes,  helicopters,  etc.,  that  form  an  excellent  background 
for  visualization.  However,  in  groundwater  flow  no  such  tangible 
objects  exist. 

c.  Heterogeneity,  in  CFD  applications  there  is  only  a  single  homoge¬ 
neous  medium  such  us  air.  Flow  in  porous  media  deals  with  signifi¬ 
cantly  different  properties  which,  at  lime.s,  varies  several  orders  of 
magnitude. 

Nevertheless,  the  following  give  a  good  representation  of  the  flow  pattern: 

a.  Color  contours.  Color  contours  on  visible  faces  used  to  show  varia¬ 
tion  of  scalar  quantities. 

b.  Isolevels.  The  3-D  equivalent  of  a  contour  plot  where  each  isolevel 
in  space  is  a  surface  representing  a  scalar  quantity  with  the  same 
value. 
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c.  Particle  traces.  Lines  in  space  that  show  the  paths  of  particles  after 
being  released  in  the  porous  medium. 

cl.  Translucency.  The  ability  to  sec  through  an  object  or  data  to  see 
what  is  behind  it. 

e.  Animation.  Viewing  many  scenes  (30  frames  per  second)  in  succes¬ 
sion  showing  movement  of  rotation,  translation,  flow,  etc. 

As  examples,  Figure  10  shows  the  hidden  line  plot  of  the  grid  with 
total  head  in  color  contours,  and  Figure  1 1  shows  the  3-D  grid  as  translu¬ 
cent  with  three  isolevels.  The  fact  that  the  isolevels  are  essentially  verti¬ 
cal  shows  that  there  is  little  variation  of  the  results  in  the  z  direction. 


Figure  10.  Total  head  contours  with  the  hidden  line  plot 
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Cray  Y-MP  improvements 


Performance  issues  vvcre  addressed  by  Mr.  Alex  Carrillo,  HPC  Center, 
WES.  For  evaluation  purposes,  in  addition  to  the  RMA  groundwater 
model,  a  smaller  aquifer  groundwater  model  and  two  different  grid  sizes 
of  the  Cerrillos  Dam  in  Puerto  Rico  seepage  model  were  u.sed  (the 
Ccrrillos  Dam  study  was  performed  by  Mr.  John  Palmerton,  GL,  WES). 
Table  1  describes  various  aspects  of  the  models.  First,  the  number  of 
nodes  and  elements  for  each  model  is  given.  Next,  the  number  of  nonlin¬ 
ear  iterations  required  to  obtain  the  steady-state  solution  is  given.  Each  it¬ 
eration  can  be  summarized  as  formulating  and  solving  the  following  set  of 
equations: 


[Kr  {Atp}'  =  -  [AQy 
{<p)'+‘  =  {(p}'  +  {A(p}' 


(3) 


Table  1 

Model  Information 

Model 

RMA 

Aquifer 

Cerrillos 

(Small) 

Cerrillos 

(Big) 

Nodes 

23,513 

11,578 

10,915 

87,572 

Elements 

38,796 

9,855 

8,469 

77,124 

Nonlinear  Iterations 

48 

5 

8 

9 

Original  Global  Bandwidth 

883 

1.594 

502 

1,974 

Original  Average  Local 
Bandwidth 

364 

1,075 

467 

1,855 

New  Globat  bandwidtti 

482 

675 

383 

1,481 

Now  Average  Local 
Bandwidth 

356 

436 

283 

1,092 

Original  CPU 

12,548 

3,747 

517 

33,140 

Original  I/O 

5,149 

910 

273 

8,040 

Now  Direct  Solver  CPU 

18a 

33 

47 

3,906 

New  Direct  Solver  I/O 

30 

1 

1 

1,235 

New  Iterative  Solver  CPU 

594 

29 

21 

278 

New  Iterative  Solver  I/O 

0 

0 

0 

0 
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where 


[K]'  -  stiffness  matrix  at  the  i'th  iteration 

(Atp)'  =  change  in  total  head  or  potential  vector  at  the  i’th  iteration 

(Ag)'  ~  residual  flow  vector  at  the  i’th  iteration 

Because  the  solution  is  steady  state,  a  significant  savings  can  be  achieved 
by  updating  [/(T]'  only  with  respect  to  changing  boundary  conditions,  and 
this  is  done.  [AQ]'  must,  however,  be  computed  using  the  current  values 
of  the  element  stiffness  matrices. 

The  initial  profiling  of  the  code  revealed  that  the  primary  area  affect¬ 
ing  performance  was  the  solution  process  (the  solution  of  the  set  of  simul¬ 
taneous  linear  equations  in  Equation  3).  Not  only  was  this  part  a 
computational  bottleneck,  but  the  I/O  required  for  the  out-of-core  solver 
also  severely  inhibited  performance.  The  I/O  for  storing  stiffness  and 
other  matrices  was  also  a  significant  factor.  Thus,  formulating  and  solv¬ 
ing  Equation  3  with  the  accompanying  I/O,  along  with  some  standard  im¬ 
provements.  were  the  primary  focus  of  the  evaluation  and  modifications. 

In  fact,  many  of  the  improvements  made  to  the  program  apply  equally 
well  to  a  generic  engineering  workstation  as  well  as  a  supercomputer. 

Solution  process.  Changes  to  the  solution  algorithm  branched  into 
two  directions.  First,  a  more  efficient  direct  solver,  a  Cholesky  factoriza¬ 
tion,  was  used  to  replace  the  Gaussian  elimination  routine  being  used  in 
the  original  program.  Secondly,  a  preconditioned  conjugate  gradient-like 
method  was  added  to  provide  an  iterative  solver  capability.  Various  prob¬ 
lem  characteristics  affected  the  performance  of  each  method.  The 
Cholesky  factorization  routine  was  developed  using  a  LAPACK  library 
subroutine  (Anderson,  et  al.  1990)  as  a  template.  Using  this  new  algo¬ 
rithm,  the  grid  was  first  divided  into  64-node  blocks.  Then,  an  out-of¬ 
core  capability  was  built  around  the  routine,  and  minor  modifications 
were  made  to  take  advantage  of  the  smaller  local  bandwidth  of  each 
block.  Also,  for  an  unchanging  stiffness  matrix  between  nonlinear  itera¬ 
tions,  reuse  of  the  factorization  from  the  previous  iteration  greatly  re¬ 
duced  the  computations  required  for  some  problems.  Coupled  with  a 
more  sophisticated  bandwidth  minimization  routine  (Gibbs,  Poole,  and 
Stockmeyer  1976),  significant  time  reductions  were  obtained  for  the  direct 
solution  process.  (Table  1  shows  botii  the  global  and  average  local 
bandwidth  using  both  the  original  and  tne  new  bandwidth  minimization 
routines.) 

The  introduction  of  the  preconditioned  conjugate  gradient-like  method 
added  an  iterative  solver  capability.  The  reduced  memory  requirements  of 
this  method  eliminated  the  need  for  an  out-of-core  solution.  Not  as  stable 
as  the  direct  solvers,  the  iterative  solvers  can  have  problems  converging 
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for  poorly  conditioned  problems.  However,  for  large,  well-conditioned 
problems,  they  can  perform  significantly  better  than  the  direct  solvers. 

I/O.  Several  changes  were  made  to  improve  the  I/O  and  memory  man¬ 
agement  of  the  program.  First,  the  stiffness  matrix  was  switched  from  a 
banded  storage  form.at  to  a  sparse  matrix  format.  This  allowed  assembly 
and  boundary  condition  modifications  to  be  completed  in-core  and  elimi¬ 
nated  many  of  the  inefficiencies  associated  with  the  out-of-core  solver. 
This  also  allowed  a  simpler  transition  between  the  iterative  and  direc.1 
solvers.  The  reduced  bandwidths  also  decreased  the  I/O  requirements 
for  the  out-of-core  solver,  making  the  use  of  faster  disks  more  feasible. 
Finally,  the  iterative  solver  eliminated  the  need  of  any  out-of-core  I/O 
time  altogether. 

General  improvements.  Enhanced  vectorization  and  other  im¬ 
provements  in  coding  were  also  accomplished.  Through  effective  use  of 
arrays,  redundant  work  was  eliminated,  and  several  functions  were  trans¬ 
formed  into  level  3  BLAS  routines  (Dongarra,  et  al.  1990). 

Results.  Table  1  shows  the  CPU  times  and  I/O  wait  times  (in  seconds) 
for  the  original  program,  the  *,ew  program  using  the  direct  (Cholesky) 
soF'er,  and  the  iterative  (preconditioned  conjugate-like)  solver.  Specific 
characteristics  of  each  problem  produced  the  performance  differences  be¬ 
tween  the  direct  and  iterative  solvers.  In  general,  however,  the  iterative 
method  is  preferred  for  large,  well-conditioned  problems.  The  RMA 
model  tended  to  produce  a  poorly  conditioned  stiffness  matrix,  resulting 
in  the  iterative  method  having  a  difficult  time  converging.  The  stiffness 
matrix  also  remained  unchanged  for  all  the  48  nonlinear  iterations,  so  sub¬ 
stantial  benefit  was  obtained  from  the  reuse  of  the  initial  factorization. 
Thus,  the  direct  solver  provided  the  best  performance.  The  aquifer  model 
was  well-conditioned,  greatly  improving  the  performance  of  the  iterative 
method.  It,  too,  benefitted  from  the  reuse  of  the  previous  factorizations 
when  using  the  direct  solver,  because  only  the  first  tv.'o  nonlinear  itera¬ 
tions  required  a  factorization.  The  net  result  was  a  comparable  time  be¬ 
tween  the  direct  and  iterative  solvers.  T  .,  Cerrillos  dam  models  were 
also  well-conditioned.  However,  in  neither  case  did  the  stiffness  matrix 
remain  unchanged,  explaining  the  superior  performance  of  the  iterative 
solver. 


2-D  Versus  3-D  Comparison 

A  very  important  aspect  of  this  research  is  to  determine  whether  a  3-D 
study  is  needed  in  the  saturated  flow  portion  of  GL’s  RMA  modeling  ef¬ 
fort.  Figure  12  is  a  color  contour  plot  showing  the  difference  in  total  head 
values  between  the  2-D  and  3-D  results.  A  3-D  result  for  a  given  2-D 
node  was  obtained  by  averaging  the  values  of  head  for  the  seven  3-D 
nodes  corresponding  to  the  2-D  node.  The  scale  is  set  this  way  so  that  the 
plot  can  be  directly  compared  to  another  plot  described  later  in  this  report. 
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Figure  12.  2-D  versus  3-D  comparison 


Yellow  is  selected  as  the  color  of  no  difference.  While  this  is  the  domi¬ 
nant  color,  the  maximum  difference  of  about  10  ft  (3.048  m)  exists  in 
some  areas.  Some  variation  is  expected,  so  this  amount  seems  reasonable. 
It  is  important  to  realize,  also,  that  the  boundary  Conditions  arc  essentially 
2-D  boundary  conditions.  Mead  is  specified  a  constant  for  the  entire  verti¬ 
cal  area  where  he;  1  is  known.  True  3-D  boundary  conditions  would  have, 
for  instance,  a  river  specified  more  accurately. 
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Grid  Display 


One  of  the  most  important  graphics  tools  to  use  in  doing  an  FEM 
model  study  is  a  program  to  display  and  edit  the  FE  grid  and  other  input 
data.  There  are  many  software  and  hardware  options  for  this,  and  the 
topic  may  even  seem  a  bit  trite  to  some.  However,  because  the  author  be¬ 
lieves  so  strongly  in  the  use  of  scientific  visualization  tools  not  only  for 
understanding  the  results  but  also  for  ensuring  the  input  data  are  correct, 
this  section  is  included.  Since  results  from  the  commercial  package 
MPGS  iave  already  been  shown,  what  is  featured  here  is  a  quickly  written 
(a  few  days)  program  called  SHOWGRID.  This  program  was  written  by 
Mr.  Scott  D.  Weberg,  SVC,  for  the  Silicon  Graphics  workstation  using 
their  Graphics  Language  software.  SHOWGRID  can  be  used  to  view  the 
2-D  grid,  its  boundary  conditions,  any  scaler  value,  and  any  negative  area 
elements.  Figure  13  shows  a  plot  of  the  menus  developed  with  initial 
head  being  plotted  for  the  node  values.  This  can  be  compared  with  the 
color  contour  plot  from  MPGS  given  in  Figure  5.  Figures  1  and  2  are  also 
SHOWGRID  plots. 

The  primary  problem  that  occurs  with  the  grid  is  negative  or  badly 
skewed  elements.  Figure  14  shows  a  typical  “bad  spot”  where  the  location 
of  an  observation  well  could  have  been  accidentally  moved  so  much  that  a 
negative  area  clement  was  created  (highlighted  element  1056).  With  pro¬ 
grams  like  SHOWGRI.O,  these  errors  can  be  easily  spotted  and  corrected. 
What  is  more  difficult,  however,  is  the  detection  of  badly  skewed  ele¬ 
ments.  Figure  15  shows  a  zoomed  area  of  the  grid  where  an  extremely 
skewed  clement  has  been  created,  and  Figure  16  shows  how  this  is  fixed. 
Numerically,  skewed  elements  such  as  these  can  cause  severe  problems, 
and  they  should  be  avoided.  However,  the  triangular  elements  turned  out 
to  be  much  less  susceptible  to  this  than  the  brick  elements  with  coincident 
nodes.  This  is  illustrated  by  the  2-D  versus  3-D  comparison  (Figure  17) 
for  a  grid  containing  several  skewed  elements  like  that  of  Figure  15.  As 
before,  one  would  expect  some  differences  in  the  two  solutions,  but  the 
rather  large  regions  where  significant  differences  occur  were  always 
caused  by  skewed  elements  in  the  3-D  grid  (the  2-D  solutions  were  essen¬ 
tially  the  same  for  the  two  cases).  In  fact,  the  large  area  where  the  differ¬ 
ence  is  as  much  as  30  ft  (9.144  m)  (dark  blue  region)  was  completely 
turned  to  yellow  by  fixing  the  one  skewed  element  shown  in  Figure  15. 
Figure  12,  corresponding  to  the  good  grid,  can  be  directly  compared  to  the 
plot  in  Figure  17  to  sec  the  dramatic  effect  of  fixing  all  such  elements. 


Summary  and  Conclusion 

The  3-D  study  showed  that  for  a  large  plan  view  flow  problem  where 
essentially  2-D  type  lH)undary  conditions  (a  constant  head  is  applied  the 
uill  depth  of  the  aquifer  and  fully  penetrating  wells  are  used)  are  em¬ 
ployed,  a  3-D  sr)Ii;tion  is  not  usually  necessary  for  the  flow  eominitations. 
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Figure  13.  SHOWGRID  menu  with  initial  head  values 


The  iO-fi  (3.048-m)  maximum  difference  in  head  is  not  significant  over 
such  a  large  area.  However,  if  one  wishes  to  apply  more  accurately  the 
boundary  condition  of,  say,  a  river  or  do  a  more  localized  study  such  as  in 
the  case  of  the  Cerrillos  Dam  analysis,  a  3-D  solution  is  required. 

Contaminant  transport  (Chapter  4)  is  another  matter.  Even  if  the  flow 
is  horizontal  and  the  heads  are  the  same  in  the  z  direction  as  computed 
from  a  2-D  plan  view  calculation,  a  3-D  study  is  required  to  get  an  accu¬ 
rate  result  for  concentration  of  contaminant.  This  is  because  advection 
will  carry  the  contaminant  horizontally  in  all  directions,  and  dispersion 
will  also  distribute  it  vertically. 
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Figure  14.  Negative  area  region 


Two-dimeiisiona!  triangular  elements  arc  more  stable  than  nine-node 
brick  elements  with  coincident  nodes  with  respect  to  skewness.  A  3-D 
prism  or  tetrahedron  element  should  be  considered,  or  a  careful  examina¬ 
tion  of  the  2-D  grid  is  required. 


This  exercise  has  also  shown  that  visualization  tools  are  essential  when 
doing  a  numerical  study  to  both  ensure  the  correctness  of  the  input  data 
and  to  properly  interpret  the  results.  Further,  without  the  use  of  a  large 
scientific  computer  such  as  the  Cray  Y-MP,  a  full  study  of  a  large  3-D 
problem  is  almost  impossible.  This  is  especially  true  of  time-dependent 
nonlinear  multiphase  flow  and  contaminant  transport  computations. 
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Figure  15.  Skewed  element 
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Figure  16.  Fixed  elements 


1 


3  Unsaturated  Flow 


Introduction 

Chapter  2  gave  some  extensive  work  using  the  FEM  for  saturated  flow. 
This  chapter  will  investigate  unsaturated  flow  using  another  computa¬ 
tional  technique  typically  used  by  the  aerospace  engineers,  which  is  the 
finite  volume  method  (FVM).  The  FVM  is  very  similar  (and  sometimes 
identical)  to  the  finite  difference  method.  The  following  will  be  presented: 

a.  Governing  equations.  The  equations  used  in  this  work  to  model  the 
unsaturated  flow  cast  in  strong  conservative  form  for  a  curvilinear 
coordinate  system. 

b.  Computational  procedure.  The  finite  volume  equations  used  to  ide¬ 
alize  the  governing  equations  and  the  numerical  scheme  used  to 
solve  them. 

c.  Comparison  of  results.  Results  obtained  from  this  FVM  technique 
arc  compared  with  FEM  results  for  various  problems,  including  labo¬ 
ratory  tests  and  analytical  solutions. 


Governing  Equations 

The  governing  partiai  differential  equation  used  for  unsaturated  flow 
for  both  the  (x,  y,  z)  and  the  (^,  t|,  Q  curvilinear  coordinate  systems  will 
now  be  given. 


(x,  y,  z)  equations 

The  equation  at  a  given  (x,  y,  z)  point  in  space  and  time  t  in  matrix 
form  is 

-I-  V^)  }  H- ^  =  F- I*  (4) 
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where 


h  =  pressure  head 

[fcj]  -  saturated  hydraulic  conductivity  matrix 
=  relative  hydraulic  conductivity 
q  =  source  or  sink 


F  =  water  capacity  given  by 


„  de 


where  6  =  is  the  moisture  content 


In  an  unsaturated  zone  of  flow,  F  and  are  functions  of  h,  making  Equa¬ 
tion  4  highly  nonlinear.  Equation  4  is  also  at  times  put  in  terms  of  the 
total  head  or  potential  <p 

<p  =  /2  +  z  -  (p^  (6 

where  <p^  =  datum 

Equation  4,  using  Equation  6,  now  becomes 


Vtp)  +  q  =  F^ 


This  work  uses  Equation  7.  Different  expres.sions  for  k^{h)  and  0(/i)  are 
used  as  these  are  empirically  based. 


(^,  q,  Q  equations 

The  approach  taken  is  to  first  map  a  region  of  geometry  into  a  (^,  r|,  Q 
curvilinear  coordinate  system.  For  example,  Figure  18  shows  three  sur¬ 
faces  of  a  Circular  region  of  soil  around  a  well  (O  type  grid  in  3-D)  that  is 
mapped  into  a  square  box  grid  in  the  (^,  q,  Q  system  as  shown  in  Fig¬ 
ure  19.  The  inner  cylinder  at  the  well  and  the  outer  cylinder  at  the  radius 
of  influence  are  q  =  constant  surfaces,  and  the  horizontal  grid  at  the  bot¬ 
tom  plotted  using  thick  lines  is  a  ^  =  constant  surface. 

Let  J  be  the  determinant  of  the  transformation  matrix  between  the  two 
coordinate  systems  for  a  nonmoving  grid  as  follows: 
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Figure  19.  Square  box  grid 

Appendix  A  contains  derivations  for  certain  equations  in  the  eurvilin- 
ear  coordinate  system.  Using  Equation  A16  for  the  strong  conservative 
form  of  the  divergence  of  some  vector  {F},  we  have 

Mif)  .  ^  (jvfiF))  +  A  ^  (jvf (fi) 

where  4,-  i  =  1. 2,  3  represent  rj,  and  C,  respectively.  Equation  9  ap¬ 
plied  to  Equation  7  produces 

3 

S  ^  j  +  Jq=JFf  (10) 


Equation  A12  can  now  be  applied  to  the  gradient  of  the  potential  to 
produce 
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Finally,  using  the  substitution 


<I>  =  y(p 


in  Equation  11,  the  following  final  result  is  obtained: 


3  3 

S  I 

«=1 ;=1 


_d_ 


dt 


Equation  1 3  is  the  equation  used  in  this  FVM  research. 


(12) 


(13) 


Computational  Procedure 

There  are  enormous  varieties  of  computational  techniques  that  can  be 
investigated  for  the  solution  of  Equation  13.  This  effort  can  only  choose 
some  of  the  most  promising  techniques  and  investigate  them.  The  chosen  so¬ 
lution  will  now  be  presented,  and  for  simplicity,  a  2-D  formulation  is  given. 


Finite  voiume  ceii 

The  standard  FV  strategy  is  to  consider  each  computational  cell  an  FV 
with  the  unknown  variables  to  be  computed  evaluated  at  the  center  of  the 
cell.  However,  the  FEM  usually  has  the  value  of  the  unknowns  computed 
at  the  grid  points  (node  points).  So,  as  shown  in  Figure  20,  a  modified 
way  of  defining  the  finite  volumes  (McCormick  1992)  is  to  use  the  dotted 


Figure  20.  2-D  computational  grid 
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lines  to  form  cells  around  the  node  points.  This  has  the  advantage  of  al¬ 
lowing  the  identical  grid  to  be  used  for  FEM  or  FVM,  but  a  disadvantage 
is  the  cells  along  the  boundary  have  reduced  size  by  one-half  or  one- 
fourth.  This  problem  has,  however,  been  minimized  by  a  judicious  soft¬ 
ware  design  of  subroutines. 


Discretized  equation 

To  understand  how  Equation  1 3  is  discretized,  first  consider  the  simpler 
equation 


3p  8<7  _  dr 

3^  3n  "  9/ 


(14) 


where  p,  q,  and  r  are  each  some  continuous  function.  Let  each  computa¬ 
tional  area  be  a  unit  square.  Also,  an  implicit  formulation  is  sought  at 
time  r  =  (n  +  l)Ar  with  an  Euler  approximation  to  the  time  partial  deriva¬ 
tive.  Then,  multiplying  Equation  14  by  d%dx\  and  integrating  over  the  unit 
square  of  a  node  at  grid  point  (/,  j)  gives 


,« +  1  +•  I  ,  n  +  1  n  +  1 


At 


(15) 


where  the  one-half  designations  represent  evaluation  of  information  on  the 
respective  edges  of  the  finite  volume  cell.  Traditionally,  it  is  stated  that 
the  value  at  the  node  point  (center  of  the  cell)  is  taken  as  the  constant 
value  over  the  entire  cell,  and  that  is  why  the  integral  signs  can  vanish. 
However,  a  better  interpretation  is  that  the  value  at  the  node  point  repre¬ 
sents  the  average  value  of  the  variable  over  the  cell.  In  like  manner, 
rather  than  thinking  of  the  one-half  terms  as  being  constant  on  the  line  seg¬ 
ments,  the  only  restraint  required  is  that  this  value  represent  the  average 
value  on  its  corresponding  line  segment. 

This  same  process  can  now  be  applied  to  a  2-D  version  of  Equation  13 
to  obtain 
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-.T 


Ji-h 


+  vn'"  ^  r/:1  f  (<PV^) 


^n  +  1  _  ^ 

+  Jq  =  I^*' 

^  '.V  At 


Let  the  remaining  partial  derivatives  in  Equation  16  be  approximated  by 
central  differences.  Then  the  first  term  of  Equation  16  becomes 


/  f  -i  1^”'*’* 
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Al'd  the  thuu  wrm  of  Equation  16  is  now 


-  Vti  «!)"■*■' 

i+i,y-i  i  +  i,_/-i 


Define  the  terms 


.  T  f  T  \  n  + 1 

{A+}T^  =  , 

I  J'  +  r-' 


{A-f  = 


[B^V  =  yn%[kU_  , 


■.j-j 


Then  Equation  16  becomes 


1  1 


>  \  c  <1»  =  -  Jq  -  —  F  <I> 

^  ^  /+m,y  +  «  j.j  Ar  i.y  /, 


1  « /I  +  I  Jtji 


m~-\  n  =  - 


where  the  C’s  are  coefficients  that  are  still  functions  of  the  dependent  vari¬ 
able,  and  thus  Hquaticn  20  remains  a  ncnlinear  equation.  NeverUieless, 
the  C’s  ill  a  nine-point  template  format  are  given  in  Table  2. 


Computation  of  geometric  quantities 

The  geometric  quantities  can  be  computed  in  a  varic:  of  ways.  Each 
element  (Figure  Al)  is  treated  as  a  four-node  isoparamcu ic  element. 
Thus,  the  transformation  between  x  and  (^,  tj)  becomes 
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Table  2 

Nine-Point  Template 

1 

2 

{B^  Tvn 

2 

2 

leilM 

2 

0 

{A'M 

-  {A*  +  A-  }'^V5 

-  {e"  +  r 

1  pl>t^ 

~  At 

-1 

LOlvii 

2 

{rfvn 

2 

^  L?  IX 

2 

n 

m 

-1 

0 

1 

x  =  (^  +  1  -  ^)  (Ho  +  1  - 11)  Xj  +  (^  -  4o)  Olo  +  1  -  n)  ^ 
—  ^)  (H  "“  ^0^  ^3  “  4)  (^  ”  ^0^  '^4 


(21) 


where 

(^.  Tlo)  =  integer  values  of  the  computational  coordinates  for 
the  first  node  of  the  computational  element 

Xi,  X2,  X3,  and  X4  =  X  coordinates  at  the  node  points 

A  similar  expression  to  Equation  21  exists  for>'.  Tenns  like 

dx 


^  =  -  (Ho  +  1  -  ^1)^1  +  (%  +  1  -  H)  ^2 
+  (H  -  Ho)^'3  -  (H  -  Ho)^4 


/">0\ 


can  now  be  computed  using  Equation  21 .  Evaluating  Equation  22  at  the 
center  of  the  computational  element  yields 


1  ,  . 
w  =  2  (''2  ■  ^ 

V 


(23) 
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Values  at  the  node  points  arc  computed  by  averaging  the  surrounding 
center  values  (four  for  internal  nodes)  using  Equation  23.  Finally,  J, 
and  Vti  are  computed  using  Equations  A4,  A6,  and  A7.  Matrices  such  as 
those  in  Equation  19  are  first  evaluated  at  the  nodes  and  then  averaged  to 
find  values  for  the  sides. 


Nonlinear  iteration 


Although  terms  in  Equation  20  can  now  be  computed,  Equation  20  is 
still  nonlinear,  so  some  iteration  scheme  must  be  adopted.  Different 
schemes  were  considered,  including  the  Picard  and  Newton  schemes  (Putti 
and  Panlconi  1992).  Presented  here  is  the  scheme  tliat  worked  well  for 
tlie  problems  tested. 


Modified  Newton  scheme.  Given  a  nonlinear  equation 

/(O)  =  0 


but  at  the  k’th  iteration 


*  0 

then  a  Taylor  Series  expansion  can  be  used  as  follows: 


=/(<!>*)  +  M  A<D  +  ... 

'k  )k 


So  a  first-order  approximation  is  to  set  the  right-hand  part  of  Equation  26 
to  zero  to  yield 


Now,  apply  Equation  27  to 


1  1 


/•(<!>*+') -y  y  .  +jq 

'  JLi  ^  i  +  jn.;  +  n  i  +  BLi  +  n  ’ 


/w  =  -l  n  =  -l 


+  —  F*'*'  ^ 
At  'w 


(from  Equation  20)  in  such  a  way  as  to  neglect  the  partial  derivatives  of  C 
and  F  to  obtain 
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(29) 


Solution  of  equations.  Equation  29  can  be  solved  by  any  number  of 
direct  or  iterative  techniques,  such  as  forms  of  Gauss  Elimination  or  the 
preconditioned  conjugate  gradient  method.  The  unknown  variables  then 
can  be  updated  v'ith  the  process  continuing  until  the  update :  become  so 
small  that  the  solution  for  the  current  time-step  has  converged.  However, 
an  under-relaxation  type  method  (the  conjugate  gradient  method  will  be 
tested  later)  worked  well  enough  for  the  problems  tested.  Because  of  the 
nonlinearity,  first  the  Gauss-Seidel  type  computation 


=  -7- 
''■>  Cf . 


1  1 


m  =  -l  n  =  -l 


+  C*  .  . 

'.y 


(30) 


was  performed.  Then,  after  results  for  all  nodes  have  been  done,  update 
using 

O*  t  *  =  ,  +  otA<l>^^  t  *  (31) 

i.y  i,j  t>j 

where  a  varies  between  zero  and  1.  Equations  30  and  31  can  be  executed 
several  times  before  updating  the  C’s  and  F,  but  it  was  more  efficient  to 
update  ail  data  after  each  iteration.  Since  only  scalar  operations  are  done, 
each  iteration  is  very  fast.  Of  course,  more  iterations  will  be  required,  so 
only  a  partial  payoff  is  realized. 


Test  Problems 


Two  test  problems  will  now  be  given  to  show  the  effectiveness  of  the 
solution  presented.  Comparison  with  FEM  results  will  also  be  made. 


Dupuit’s  problem 

The  first  problem  is  the  classic  problem  of  steady  state,  unconfined 
flow  in  an  earth  embankment  with  vertical  sides  and  an  impervious  base 
as  shown  in  Figure  21.  Water  flows  from  headwater  to  the  exit  face  where 
the  pressure  =  0  line  represents  the  free  surface.  The  embankment  is  100 
by  100  ft  (30.48  m),  the  headwater  level  is  at  100  ft,  and  the  tail  water  is  at 
20  ft  (6.096  m).  The  exit  point  where  the  free  surface  intersects  the  down¬ 
stream  vertical  boundary  is  40  ft  (12.192  m).  A  grid  of  1 1  X  1 1  =  1 21  nodes 
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Figure  21.  Dupuit’s  problem 

and  10  X  10  =  100  elements  with  Ax  =  10  ft  (3.048  m)  and  Ay  -  10  ft 
(3.048  m)  was  built  and  run  on  the  developed  FV  program  and  a  2-D 
saturated,  confined/unconfined  flow  FEM  program.  Again,  because  of  the 
way  the  FV  program  was  constructed,  the  identical  grid  and  boundary  con¬ 
ditions  can  be  used  in  both  codes. 

The  2-D  FEM  program  (2DSEEP)  computes  “stiffness”  type  matrices 
and  then  assembles  them  to  form  the  set  of  equations  to  be  solved.  These 
stiffness  matrices  are  numerically  integrated  by  evaluating  the  integrand 
at  certain  integration  points  and  then  using  a  Gauss  quadrature  formula. 
Typically,  2x2  =  4  integration  points  are  used,  but  because  the  free  sur¬ 
face  can  go  through  an  element,  2DSEEP  uses  4x4=  16  integration 
points.  The  relative  hydraulic  conductivity  is  then  set  to  0.(X)1  for  an  inte¬ 
gration  point  with  pressure  head  of  less  than  zero.  This  integration  pro¬ 
vides  a  smearing  process  so  the  transition  is  not  so  abrupt. 

Since  no  integration  is  used  in  the  FV  algorithm,  a  5-ft  (1.524-m)  tran¬ 
sition  zone  was  therefore  provided  where  the  relative  hydraulic  conductiv¬ 
ity  varied  from  1  to  0.001  when  the  pressure  head  varied  from  0  to  -5  ft 
(-1.524  m).  Relative  hydraulic  conductivity  was  then  kept  at  0.001  for  all 
values  of  pressure  head  less  than  -5  tt  (-1.524  m). 

Finally,  the  moisture  capacity  F  in  Equation  5  is  set  to  zero  to  eliminate 
the  time-dependent  aspect  of  Equation  4  for  a  steady-state  solution. 

Figure  22  shows  a  comparison  of  the  free  surface  for  the  FEM  and  FV 
solutions.  Note  that  they  are  rather  close.  Also,  the  running  times  of  the 
two  solutions  can  be  compared.  It  must  be  emphasized  that  the  FEM 
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Free  Surface 


Figure  22.  Comparison  of  FEM  and  FV  free  surface 


technique  has  significantly  more  computations  I  jcause  of  the  numerical 
integration  process  and  the  solution  of  a  banded  system  of  equations  at 
each  iteration.  In  general,  fewer  elements  for  the  FEM  .solution  are  there¬ 
fore  needed  compared  to  the  FV  grid.  Also,  because  a  banded  system  is 
solved  compared  to  a  relaxation  iteration,  much  fewer  FEM  iterations  are 
needed.  With  these  thoughts  in  mind,  the  FEM  solution  took  1  min, 

1 3  sec,  and  seven  iterations  on  a  486  class  PC  running  at  33  mhz,  and  the 
FV  solution  took  12  sec  with  124  iterations.  It  should  also  be  observed 
that  the  diagonal  terms  in  the  template  in  Table  2  go  to  zero  for  this  square 
grid,  making  the  formulation  similar  to  a  standard  finite  difference  algorithm. 


Laboratory  test  problem 
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Results  from  an  experimen  il  study  of  2-D  transient  un.saturated/saturated 
flow  with  water  table  recharge  (Vauclin,  Khanji,  and  Vachaud  1979)  will 
now  be  compared  with  results  from  both  the  FV  formulation  and  a  3-D 
transient  unsaturated/saturated  FEM  flow  code.  The  problem,  as  shown  in 
Figure  23,  consists  of  flow  i.i  a  homogeneous  soil  of  saturated  hydraulic 
conductivity  35  cm/hr  in  a  tank  600  cm  long,  200  cm  tall,  and  5  cm  thick 
with  an  impervious  bottom.  Because  of  symmetry,  only  300  cm  of  tlie 
tank  are  modeled  with  the  center  line  (AF  in  Figure  23)  being  treated  as 
an  impervious  boundary.  A  constant  pool  elevation  of  65  cm  is  main¬ 
tained  along  BC  in  Figure  23  with  the  boundary  CDE  covered  to  avoid 
evaporation,  EF  i.s  initially  covered  and  the  tank  is  allowed  to  completely 
settle.  Then  EF  is  uncovered  and  a  flow  rale  of  14.8  cm/hr  is  applied  to 
the  system  for  8  hr  while  holding  BC  at  the  constant  total  head  of  65  cm. 
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Figure  23.  Laboratory  problem 

The  relative  hydraulic  conductivity  in  the  unsuturated  n':.;!  .n  was  deter¬ 
mined  experimentally  to  be 


/I  +  (-/■) 


where 


-4  =  2.99  X  10^ 

B  =  .‘i.O 

The  moisture  content  equation  was  also  determined  experimentally  and  by 

t  '1  I  «• t  l>«a 

VIIV  UUV  \Jt  U  IW44i>k~  Sr>J  l-vr 


0  =  0 


a  +  (-/i)P 


where 


0,  =  0.30 


a  =  40.00 


p  =  2.90 
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The  grid  (shown  in  Figure  24)  consists  of  a  16  x  16  structured  grid 
with  the  intervals  slightly  different  to  align  with  key  points.  For  example, 
the  dot  is  a  node  at  (161,  79)  where  pressure  head  data  were  collected.  At 
was  set  to  0.05  hr  and  allowed  to  grow  20  percent  a  time-step  until  a  maxi¬ 
mum  of  At  =  1  hr  was  reached.  Twenty  time-steps  were  done  for  a  total  of 
8  hr.  The  proper  subroutine  was  modified  to  incorporate  Equations  32, 

33,  and  the  derivative  of  Equation  33  (from  Equation  5),  and  the  problem 
was  run  on  both  the  FEM  and  FV  codes.  Because  of  the  nature  of  the 
water  capacity  curve  (Figure  25),  the  FEM  program  would  not  converge, 
and  the  FV  solution  converged  only  by  using  a  small  a  =  0.1  in  Equa¬ 
tion  31  for  the  first  few  time-steps.  The  FEM  algorithm  uses  a  Picard- 
type  iteration  strategy,  while  the  FV  scheme  uses  the  Newton-type 
iteration  as  previously  described.  However,  tabular  forms  of  Equa¬ 
tions  32,  33,  and  F  (27  data  points),  with  F  not  allowed  smaller  than  0.001 
for  large  -h,  were  used  with  the  FEM  code,  and  convergence  was  then 
achieved.  A  Picard-type  algorithm  was  implemented  in  the  FV  code,  and 
the  same  lack  of  convergence  was  observed.  Newton-vcrsus-Picard  itera¬ 
tion  is  an  important  research  topic,  so  at  least  for  this  problem  and  im¬ 
plementation,  the  Newton  algorithm  is  superior. 


Figure  24.  Grid 

The  results  for  the  40  time-.steps  for  both  the  FV  and  FEM  solutions 
were  then  obtained,  and  the  free  surface  was  compared  with  the  laboratory 
results  as  shown  in  Figure  26.  The  dis.sipativc  error  in  the  FEM  solution 
is  more  than  that  of  the  FV  implementation.  F  being  modeled  more  accu¬ 
rately  in  the  FV  code  could  have  had  an  impact  as  well. 
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Water  Capacity 


Figure  25.  Water  capacity  curve 


Figure  26.  Comparison  of  laboratory,  FEM,  and  FV  results 
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Summary  and  Conclusions 


There  are  many  more  problems  that  can  be  considered  in  unsaturated 
flow.  There  are  problems  with  difficult  geometries  and  heterogeneities. 

In  some  cases,  the  FEM  will  function  better,  and  in  other  cases,  as  has 
been  demonstrated,  the  FVM  will  function  better.  But  it  can  be  concluded 
on  the  basis  of  this  study  that  alternative  techniques  such  as  the  FVM  can 
be  a  powerful  tool  for  groundwater  modeling. 
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Introduction 

Chapter  3  showed  the  application  of  the  FVM  to  unsaturated  flow  with 
good  results.  This  chapter  will  investigate  the  ef  fectiveness  of  the  FVM 
for  contaminant  transport.  There  are  several  higher-order  numerical  tech¬ 
niques  with  one  being  the  Lagrangiau-Eulerian  Method  (LEM)  (Yeh 
1990).  This  technique  has  almost  exclusively  been  used  with  the  FEM,  so 
this  chapter  shows  the  results  of  R&D  using  both  FVM  and  LEM.  The 
following  will  be  presented: 

a.  Governing  equations.  The  equations  used  in  this  work  to  model  con¬ 
taminant  transport,  including  the  one  used  for  the  LE  approach. 

b.  Computational  procedure.  The  LE  algorithm  in  an  FV  environment. 

c.  Comparison  of  results.  Results  obtained  from  this  FVM  technique 
are  compared  with  F.HM  results  for  analytical  solutions. 


Governing  Equations 

The  governing  partial  differential  equation  for  contaminant  transport  in 
unsaturated  porous  media  at  a  given  (x,  y,  z)  point  in  space  and  time  r  is 

+  pftf  +  V  .  (vq  =  V  .  (0Z)  .  vq 

-  X(QC  +  p^5)  +  QC.^ 

where 

0  =  moisture  content 

C  =  material  concentration  in  aqueous  phase 
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=  bulk  density  of  the  medium 
S  =  material  concentration  in  adsorbed  phase  {MIM) 

V  =  discharge  velocity  vector 
D  =  dispersion  coefficient  tensor 
X.  =  decay  constant 
Q  -  source  rate  of  water 
C,„  =  material  concentration  in  the  source 
For  this  study,  the  linear  isotherm  model  for  adsorption  is  used,  which 


=  Kf 


where  Kj  is  the  distribution  coefficient 


The  component  of  the  dispersion  coefficient  tensor  d^j  is  given  by 


=  e 


“r  -  “i-* 


V.V. 

-U 

ivl 


+  a  t6.. 

m  ij 


is 

(35) 


(36) 


where 


Uf  =  lateral  dispersivity 
ai  =  longitudinal  dispersivity 
V,-  =  j'th  component  of  v 
a„,  =  molecular  diffusion  coefficient 
T  =  tortuosity 


6.J  =1  i=j 


(37) 


=  0  i 

Equation  34  can  be  expanded  and  combined  with  Equation  35  to  produce 


dC 


dt 


+  V  •  V 


(38) 


=  V  •  (0D  •  VC)  -  X{Q  +  C  +  QC. 


44 


Chapter  4  Contaminant  Transport 


Also,  unsaturated  flow  is  modeled  by 


V  ■ 


V  =  Q 


(39) 


Equations  38  and  39  can  be  combined  to  produce  the  following  governing 
equation  used  in  this  study; 


j'e  +  ~  +  V .  VC  =  V  •  (ez>  •  VQ 

-[X('0  +  p,^:^j+S]c+  ec, 


Euierian  Versus  Lagrangian  Approach 


The  typical  approach  in  2-D  is  to  stay  at  a  given  (jc,  y)  position  in  space 
and  observe  phenomena  as  they  pass.  This  is  known  as  the  Euierian  sys¬ 
tem  and  is  represented  by  Equation  40.  Another  approach  is  to  have  the 
observation  point  move  as  well  at  the  same  velocity  as  the  particles  of 
fluid  (Lagrangian  system).  The  Euierian  approach  is  modeled  using  a 
fixed  grid,  and  the  Lagrangian  approach  can  be  handled  using  a  moving 
grid  with  each  node  point  moving  at  their  respective  velocities.  This  can 
lead  to  a  skewed  grid  after  a  while,  so  a  modified  approach  (LEM)  is  used 
in  this  study. 

Using  the  chain  rule  of  differentiation,  the  change  of  concentration  for 
a  node  in  a  moving  grid  is  given  by 


A)  + 
X 


(41) 


where  the  subscript  sp  represents  grid  point,  the  subscript  fp  represents  a 
fixed  (x,  y)  point  in  space,  is  the  speed  of  the  grid  point  in  the  x  direc¬ 
tion,  and  is  the  speed  of  the  grid  point  in  the  y  direction.  The  first  time 
derivative  is  often  written  with  a  capital  D,  and  the  second  time  derivative 
is  what  appears  in  Equation  40,  so  rewrite  Equation  41  as 


DC 

Dt 


3^  +  o  ■  VC 
3/ 


(42) 


where  o  =  grid  point  velocity  vector 

Equation  42  can  now  be  substituted  into  Equation  40  to  get 
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45 


(9  +  PA) 


DC  ^ 
Dt  ^ 

\ 


e  + 


-  -u  •  VC 


(43) 


=  V  -  (az)  •  VO  -  +  e] c  +  ec. 

So  if  the  velocity  of  the  grid  point  is 


u  = 


(44) 


0  + 

then  the  Lagrangian  equation  representing  contaminant  transport  becomes 
DC 


0  + 


(.QD.vq 


(45) 


-  rxfe  + 

L 


Pc*^^)  +  e] 


C  +  QC. 

^  in 


Computational  Procedure 


It  must  be  emphasized  that  the  resulting  concentrations  from  Equation 
45  are  for  a  grid  that  has  moved.  Therefore,  the  following  two  steps  must 
be  performed  in  either  order; 


a.  Lagrangian.  Compute  an  advection  contribution  using  grid  point 
speeds  according  to  Equation  44. 

b.  Eulerian.  Compute  dispersion  and  other  contributions  using  Equa¬ 
tion  45. 


In  this  work  the  advection  step  is  first. 


Lagrangian  step 

The  FV  transport  code  that  was  written  for  this  research  takes  as  input 
the  moisture  content  and  discharge  velocity  for  each  node  and  each  time- 
step.  Therefore,  for  a  given  node  P  (as  shown  in  Figure  27  where  both  the 
velocity  components  are  positive),  grid  speed  velocities  at  time-step  n  can 
be  computed  by 

„  (46) 
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Figure  27.  Grid  point 

The  distance  node  point  P  travels  during  a  time-step  Ar  is 
Ar  =  ^  -f-  u"  '-']  At 


(47) 


Now  if  the  approximation  is  made  that  point  P'  travels  the  same  Ar  to 
reach  point  P  for  the  time  interval  Ar  (backward  tracking),  then  the  advec- 
tion  concentration  C  at  time  ri+1  for  node  P  is  simply  the  concentration 
at  P'.  However,  the  grid  velocity  at  P'  is  not  in  general  the  same  as  it  is  at 
P,  leading  to  possible  error.  A  correcting  formulation  follows. 


Computation  of  C*  at  P'.  As  before,  assume  a  four-point  isoparamet¬ 
ric  finite  element  formulation  for  each  cell  as  shown  in  Figure  A1 ,  except 
that  now  let  the  parameters  T]  vary  between  0  and  1.  Then  inside  each 
element 


X 

y 

c 

■  =  (1  -  ^)  (1  -  Tl)  • 

X 

y 

c 

>  +  ^(1  -  n)  • 

X 

y 

c 

o 

['>J 

1 

V 

(48) 


+  ^ri 

X 

y 

c 

+  (1 

-  ^)T!  • 

X 

y 

c 

V 

3 

M 

The  more  accurate  equation  for  /,  the  position  of  P',  using  Equation  47  is 


2 

2 


+  •u'''  +  ’jAr  = 


(49) 
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where 


\)'  =  grid  speed  at  P' 

rp  =  vector  containing  the  (jr,  y)  coordinates  of  point  P 
For  the  <’th  node  of  the  cell,  define  the  quantities 


Applying  Equation  48  to  the  prime  variables  in  Equation  49,  using  the  def¬ 
inition  of  Equation  50,  and  collecting  terms  yields  a  set  of  equations  to  cal¬ 
culate  Ti')  as  follows: 

(a:,  -  4-  ;^3  -  X^)  ^'Ti'  +  (X^  -  xp 


+  (X4  -  X,)ti'  =  A—  X, 


(Kj  -y^  +  Y^-  Y^)^'^'  +  (Y^  -  F,)4' 


2  '  1^ 


Rewrite  these  equations  using  constants  etj,  t/j*  ^2’  ^1’  ^2'  ^2 

a,4'Ti'  + 

(52: 

-.'•n'  +  +  CjTj'  =  ^2 

1!,'  first  of  these  equations  can  be  solved  for  Ti'  to  yield 

-  b^^:  (53; 

’i'  =  — 1 — i' 

c'l  + 

Equation  53  can  now  be  substituted  into  the  bottom  portion  of  Equa¬ 
tion  52  to  obtain  the  quadratic  equation 

{a^b^  -  -  {a^d^  -f  b^c^  -  c^b.^ 


-  dM^%'  -  c,d^  +  d.c,  =  0 


which  can  now  be  easily  solved  for^'.  This  done,  Tj'  can  be  computed 
from  Equation  53,  and  finally,  substituting  the  values  of  and  q'  into  the 
C  row  of  Equation  48,  the  value  of  concentration  at  P'  is  given. 
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Determining  what  cell  to  locate  P\  Figure  27  shows  that  when  both 
components  of  grid  speed  are  positive,  the  negative  x  and  y  direction  is 
where  the  proper  cell  is.  Since  this  FV  formulation  is  on  a  structured  grid, 
the  process  of  finding  the  proper  cell  is  to  simply  look  in  the  four  adjacent 
cells  by  solving  Equation  51  for  and  T)'  and  selecting  the  cell  where 
and  Tj'  arc  both  between  0  and  1. 


Eulerian  step 

The  final  step  is  to  use  the  C*  values  in  conjunction  with  an  FV  repre¬ 
sentation  of  Equation  45  to  determine  the  concentration  at  each  node. 
Since  Equation  45  is  remarkably  similar  to  Equation  7,  a  very  similar  for¬ 
mulation  can  be  used.  With  the  definition 


o  =  yc 

the  curvilinear  coordinate  version  of  Equation  45  is 


(55) 


3  3 


S  I 


-  iX(e  +  +  Q\  <j> 


(56) 


4-  JQC„  =  (0  4. 

The  I'V  representation  i.s 


H-  I 


/  ...f  0 


7'  +  2.y 


at 


1  , 


1 T  ‘ 

U'-L' 


H- 


ViV.{0|z;i^^(<i)V^) 


yi  1  ^ 

.  ,  1 
J>-J+  0 


Vri' joiD]  1^-  ((I)V4)[ 


Y<  + 1 


(57) 


^  [  7  1  Y‘  '  ' 

Vr/ |0[z;]  (d>Vii)|  J 


/ 


a 


1  v  + 


Vii' ^0[D) 


A7- 


.  1 


<1)"  '  -  <I>* . 

-  f  Xa"  '  4-  (2  W"  ‘  4-  JQC.  =  a'.'  ^  - ^ 

1  '.V  ^  III  i-j  A/ 
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where 


a'lj  *  =  Q'I'J  ’  +  (58) 

It  is  important  to  note  that  4»*  appears  in  the  time  derivative  term  instead 
of  <I)".  The  exact  same  solution  technique  used  previously  can  be  used 
here,  except  now  the  system  of  equations  is  linear. 


Test  Problems 

Two  test  analytical  problems  will  now  be  given  to  show  the  effective¬ 
ness  of  the  solution  presented.  Comparison  will  be  made  with  FEM 
results. 


1-D  steady-state  transport 

The  problem  consists  of  one-dimensional  (1-D)  steady-state  contami¬ 
nant  transport  over  a  length  L  in  a  flow  of  discharge  velocity  u,  a  moisture 
content  e,  and  a  dispersivity  D.  The  initial  concentration  C  of  the  contami¬ 
nant  is  0,  and  then  a  spill  occurs,  giving  boundary  conditions  of  C  =  1  at 
a:  -  0  and  C  =  0  at  jr  =  L.  The  governing  equation  is 


dx^  0  dx 


(59) 


The  value  of  the  variables  used  are  L=  100,  «  ==  0.2,  0  =  0.4,  D  =  50, 
and  A/  =  10.  A  grid,  11x2,  was  created  with  a  constant  Ajc  =  10,  and  the 
problem  was  run  on  both  the  FV  code  and  an  LE  FEM  program.  The  FEM 
and  FV  results  were  essentially  identical,  and  a  plot  of  the  analytical  ver¬ 
sus  numerical  results  is  given  in  Figure  28.  These  results  are  certainly  ac¬ 
ceptable.  However,  if  u  is  increased,  a  finer  or  adaptive  mesh  is  required 
for  the  same  quality  result. 


2-D  transient  transport 

This  problem  consists  of  saturated  flow  in  a  rectangular  vertical  cross 
section  of  sand  of  size 
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Figure  28.  1-D  steady-state  results 


0  ^  X  <  a 
0  ^  y  <  L 


(61) 


that  is  initially  clean  until  a  spill  occurs  on  the  top  of  the  sand  (y  =  L).  A 
concentration  Cq  in  a  strip  of  length  s  in  the  middle  of  the  top  of  the  sand 
is  maintained  for  a  time  Iq,  and  then  it  decays  expedientially  with  a  decay 
constant  a.  Water  is  flowing  in  the  -hot  direction  with  a  discharge  velocity 
u.  However,  no  contaminant  due  to  dispersion  flows  through  the  bound¬ 
ary  aix  =  a.  Adsorption  into  the  medium  of  bulk  density  occurs  lin¬ 
early  with  a  distribution  coefficient  of  K^. 


The  governing  partial  differential  equation  is 


+  =  e  [d, 

dx  s  1 


(62) 


where 


=  maximum  moisture  content 
D\  =  dispersion  coefficient  in  the  x  direction 
D'l  =  dispersion  coefficient  in  the  y  direction 
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The  initial  conditions  are 


Chaptoi  ■  mtamlnani  Transport 


Also, 


d.  = 


f%  «  A 

a  - 

p  sink^ 

2 

1 

. 

V  y 

+  cosX^ 


—  i'' 


2 

V  y 


-  e 


a  -  s 

-P  2 


P  sinX. 


— 2— 


'  2 

V  y 


The  two  transient  terms  7’|  and  T2  are  given  by 


(69) 


T  =  -■-  (1  _  ,-P') 

‘  JA 


(70) 


and 


^2  ~ 


11  a  -  It  11(11  -  <x) 


U(t  -  /(,) 


(71) 


whore 


R 


2  • 

D,  r 4-  Xf]  4- 

'nn 

~L 

-  V  J 

K  J 

•J 

and 


U(t  -  to)  -  '  <  \) 

=  1  1  ^  r„ 


(72) 


(73) 


The  value  of  the  parameters  are  given  in  Table  3.  A  21  by  1 1  grid  was 
used  to  solve  the  problem  using  both  the  TliM  and  bV  codes,  bigure  29 
shows  a  comparison  of  analytical,  bBM,  and  b'V  results  for  position  (100, 
1 8)  for  10  time-steps  where  At  -  1.0.  A  higher  level  of  accuracy  can  be 


Tabic  3 

Value  of  Parametets 

Co 

1.0 

a 

100.0 

L 

20.0 

s 

20.0 

to 

10.0 

a 

1.0 

0, 

0.4 

po 

1.2 

Ko 

0.1 

u 

2.0 

Ui 

60.0 

Di- 

5.0 
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Figure  29.  2-D  transient  results 

uchiuvcd  by  refining  the  grid  near  the  tup  where  the  concentration  is  ap¬ 
plied.  However,  the  current  grid  is  sufficient  to  show  the  desired  compari¬ 
son  of  results. 


Scientific  Visuaiization 

Visualization  of  contaiuinant  transport  analysis  results  is  an  important 
aspect  of  understanding  the  volumes  of  output  data.  There  arc  .several 
techniques  used,  and  there  is  potential  for  more  research,  especially  with 
the  use  of  transluccney.  One  technique  is  to  take  a  value  of  the  contami¬ 
nant  and  creatf  the  3  D  geometry  representing  the  isolcvel  and  perform  a 
hidden  surface  plot.  Another  option  is  to  cut  a  cross  .section  through  the 
plume  at  a  particular  time  and  plot  a  line  or  color  contour  of  concentration 
on  the  cross  section.  Figure  30  shows  a  plume  from  a  3-D  version  of  the 
above  problem  at  r  =  10.  The  spill  now  occurs  on  an  s  by  s  square  arcu  on 
the  top  of  the  u  X  bx  L  rectangular  region  of  porous  media  with,  in  this 
ease,  a  =  h  ^  100.  F'iguie  30  contains  an  isolevel  of  C  =  0.4,  a  color  con¬ 
tour  plot  of  concentration  along  the  vertical  cross  .section,  x  =  50,  and  a 
color  contour  plot  of  concentration  along  the  vertical  cross  section,  y  =  50. 

As  shown  in  Figure  1 1  for  saturated  flow,  the  surrounding  porous  media 
can  be  represented  using  translucent  colors, 
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|-igiiio  30.  li'.olovol  and  coloi  conlouiT. 

Summary  and  Conclusions 

As  in  lln'  l  asa  of  iin.s,iluiali.'d  I'luw,  lliiTi'  aiv  iiiiiiiy  niuiv  probk'iii.s  dial 
can  In.'  cnii.siilcicil.  'I'liciv  arc  piuldciii.s  widi  ililTiciill  j'cnniclrics  and  hclcr 
oj’cnciiics.  In  sonic  cases,  die  l  liM  will  rinicdoii  heller,  aiul  in  oilier 
eases,  as  lias  been  ilenioiislraled,  die  b'VM  will  liinelion  heller.  Hiil  il  can 
be  r'oneliided.  as  belore,  on  the  li.isis  ol  iliis  sliuly  dial  alleriiali'ii’  Uwli- 
nn|iie.s  such  as  ilu'  I'VM  can  be  a  powcn'iil  lool  lor  ’'rouiulwaU'i'  nUHleliii;;. 


(ai.ipli’ra  Conl.iniiii.inl  Ii.ims-|)o:1 
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Appendix  A 

Derivation  of  Curviiinear 
Coordinate  Equations 


Purpose 

The  purpose  of  this  appendix  is  to  derive  expressions  used  to  east  par¬ 
tial  differential  equations  in  strong  conservation  form  for  a  curvilinear  co 
ordinate  system.  These  include  the  geometric  conservation  law,  the 
gradient  operator,  and  the  divergence  operator.  2-D  versions  of  the  equa¬ 
tions  are  derived  for  simplicity. 


Coordinate  System 

Each  quadrilateral  cell  or  element  of  a  finite  volume  or  l  inite  element 
grid  is  transformed  into  a  square  (see  Figure  A1 )  fioiii  an  ( v)  coorilinat! 
system  to  a  curvilinear  (^,  t])  coordinate  system. 


3 


2 


Figure  A1.  Transformation 
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Geometric  Conservation  Law 


The  geometric  conservation  law  is  a  basic  tool  for  placing  equations  in 
strong  conservative  form  (Thompson,  Warsi,  and  Mastin  1985).  A  matrix 
oriented  derivation  of  this  equation  will  now  be  given.  The  chain  rule  of 
differentiation  states 

if  _  .  ifiz 

dx  dy 

where /is  a  continuous  function.  The  matrix  version  of  this  equation  is 

a/  iz  if  (A2) 

.a^  L  a^^. 

if  “  a^  ^  _a/ 

,ariJ  [an  anj  [a>-_ 

An  alternate  way  is 

a/l  ^  ^  fi/^l  (A3) 

ajc  _  ajc  ajc  a^ 

'if r  ii  ^ 

.ayj  [ay  ayj  [an. 

Let  the  determinant  of  the  matrix  in  Equation  A2  be  defined  by 

.  ii.ay  (A4) 

^  a^  dn  ■  an  a^ 

The  inverse  of  the  matrix  in  Equation  A2  can  also  be  multiplied  to  both 
sides  of  Equation  A2  to  produce 


ijL  _  iz]  [if 

an  a^  a^ 

dx  dx  df 

an  a^  [an 


Comparing  Equations  A3  and  A5  gives 

'i^l  iz' 

IrH  I' 

a>.  l  an. 

*  Refeiciiccs  cited  in  this  appendix  arc  listed  following  the  inuin  icxt. 
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and 


dx 

1 

>  =  —  ^ 

j 

dx 

ad 

Equations  A6  and  A7  can  now  be  used  to  evaluate 


a(W^)  acyvq) 

.  ^  . 

3^y 

“  a^an 

85  8rr  ■' 

d^x 

^x 

a^.anj 

a^anj 

It  is  clear  now  that 


nm  3(jvii) 

85  +  8n  "" 

which  is  the  geometric  conservation  law. 


(A7) 


(A8) 


(A9) 


Gradient 


The  gradient  of/from  Equation  A3  is 


(A  10) 


Multiplying  Equation  A 10  by  7  and  Equation  A9  by  / and  adding  gives 


(All) 


Combining  terms  gives  the  desired  result 


(A  12) 


Divergence 


The  divergence  of  a  vector 


(/■-}  = 


(A13) 
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A3 


in  matrix  notation  becomes 


V^F)  = 


3/.  3/, 


(A  14) 


dx  dy 

Now  if  Equation  A12  is  modified  to  the  form  of  an  operator,  this  yields 


yV  o  =  ^  (JV^  o)  +  -^  (JVrt  ° ) 

Taking  the  transpose  of  Equation  A15  and  operating  on  (F)  gives 
yv^{F}  =  (yVTi^{F) ) 

which  is  the  desired  result  for  the  divergence. 


(A  15) 


(A  16) 


A4 


Appendix  A  Derivation  of  Curviiinear  Equations 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

0MB  No  0704-0188 

Tublic  reporting  burden  for  thi»  ccHiKtion  of  inform«uon  ii  estimated  to  average  \  Kour  per  response.  fn<luding  the  time  for  reviewing  instructions,  searching  existing  data  sources. 
gat^erlng  and  maintaining  the  data  needed,  and  completing  and  reviewirvg  the  coUection  of  information  Send  comments  regarding  this  burden  estimate  or  any  other  afpect  of  this 
collection  of  information,  including  suggestions  for  r^ucing  this  burden,  to  Washirsgton  Headquarters  Services,  Directorate  for  information  Operations  and  Repoas,  Ul5  Jeffecson 
Davis  Highway.  Suite  1204.  Arlington.  VA  22202*4)0].  arvj  to  the  Office  of  Management  and  Oc^gct.  Paperwork  Reduction  Project  <0704*0  t|$).Washirsgton,  DC  20S03. 

1.  AGENCY  USE  ONLY  (L««v«  blink)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  OATES  COVERED 

September  1994  Final  report 

*  ^XppftcatKm  oPt^omputational  and  Visualization  Methods  to 

Groundwater  Modeling 

5.  FUNDING  NUMBERS 

t.  AUTHOR(S} 

Fred  T.  Tracy 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORESS(£S) 

U.S.  Army  Engineer  Waterways  Experiment  Station 

3909  Halls  Ferry  Road,  Vicksburg,  MS  39180-6199 

B.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

Technical  Report  ITL-94-7 

9.  SPONSORING /MONITORING  AGENCY  NAM£(S)  AND  ADORESS(ES} 

Assistant  Secretary  of  the  Army  (R&D) 

Washington,  DC  20315 

10.  SPONSORING /MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

Available  from  National  Technical  Information  Service.  5285  Port  Royal  Road,  Springfield,  VA  22161. 

IZa.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

_ 

13.  ABSTRACT  (Maximum  200  words) 

This  report  documents  basic  and  applied  research  in  tlic  area  of  groundwater  modeling.  A  three-dimensional 
finite  clement  method  program  and  model  for  the  real-world  problem  of  flow  near  Rocky  Mountain  Arsenal 
(RMA)  is  described,  and  a  comparison  with  a  two-dimensional  finite  clement  method  RMA  compulation  is 
done  to  help  determine  when  a  three-dimensional  analysis  is  needed.  Scientific  visualization  techniques  and 

Cray  Y-MP  improvements  developed  in  this  RMA  study  are  described.  The  finite  volume  method  is  applied 
to  saturatcd/unsaturaled  flow  and  contaminant  transport,  and  the  results  arc  compared  to  finite  clement  method 
solutions,  Nonlinear  Newton  iteration  for  unsatuiatcd  How  and  tlte  Lagrangian-Eulcrian  method  for  contami¬ 
nant  transport  arc  discussed.  Experimental  results  and  analytical  solutions  are  used  to  compare  different  proce¬ 
dures.  Finally,  contaminant  transport  results  arc  viewed  using  scientific  visualization  tools. 

14.  SUBJECT  TERMS 

Containment  transport  Groundwater  modeling 

Finite  element  nietliod  Scientific  visualization 

Finite  volume  method  Unsaturated  flow 

15.  NUMBER  OF  PACES 

68 

ie.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

UNCLASSIFIED 

18.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

UNCLASSIFIED 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRAO 

20.  LIMITATION  OF  ABSTRACT 

NSN  7540-01-280-5500 


Standard  Form  298  (Rev  2-89) 

Pr^rib«d  bv  ANSI  Std  £J9-l& 

^9&-102 


